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We study causal viscous hydrodynamics in the context of central relativistic heavy-ion 
collisions and provide details of a straightforward numerical algorithm to solve the hydrody- 
namic equations. It is shown that correlation functions of fluctuations provide stringent test 
cases for any such numerical algorithm. Passing these tests, we study the effects of viscosity 
on the temperature profile in central heavy-ion collisions. Also, we find that it is possible 
to counter-act the effects of viscosity to some extent by re-adjusting the initial conditions. 
However, viscous corrections are strongest for high-mass particles, signaling the breakdown 
of hydrodynamic descriptions for large n/s. 



I. INTRODUCTION 

Successful fits of ideal hydrodynamics to experimental data on several observables [l|, 2, 0, 0, 0] 
at the highest energies of the ongoing heavy-ion program 

aaaa 

at the Relativistic Heavy-Ion 

Collider (RHIC) seem to indicate a very small value of the ratio shear viscosity over entropy. Since 



calculations of this ratio in QCD at weak coupling a s <C 1 give Hj, lll( a numerical value that 



turns out to be larger by about one order of magnitude than the conjectured strong coupling value 




), this has given rise 

3, MM- 



for relativistic quantum field theories at finite temperature [12 [_ 
to the idea of a "strongly-coupled" quark-gluon plasma phase 

However, up to now calculations in the regime of strong coupling are limited to theories which 
possess a gravity dual theory, such as M = 4 Super Yang-Mills Theory [2l|. So far, no such dual 
theory has been discovered for QCD, meaning that the main available tool to study dynamical 
processes in QCD are based on weak-coupling approaches (although lattice-based techniques for 
making quantitative measurements of near-equilibrium quantities may be available soon 22. |23|). 

Given that the numerical value of the QCD coupling a s within the range of temperatures 
applicable for RHIC is assumed to be close to the range a s = 0.2 — 0.4, one observes that while 
this value is not very small, it is not very large either. Thus, although it would be of great 
interest to have results for QCD at very strong coupling, there is at least some hope that existing 
weak-coupling techniques might actually offer a description of RHIC physics that is not inferior to 
still-to-be-discovered QCD strong-coupling techniques (or likewise, extrapolating existing strong- 
coupling results from theories (very) different than QCD). 

Along these lines, it has recently been discovered that non-Abelian plasma instabilities [24| 
create turbulent color magnetic fields 25j, [26|] that may induce a very small effective ( "anomalous" ) 
shear viscosity coefficient [27I . |28| . without invoking strong coupling effects. Within the initial 
conditions obtained in the color-glass-condensate model [29] it is, however, unclear whether this 
effect is relevant for present RHIC energies 3y, |3l|] . 

Regardless of these issues, it is important to note that so far the ratio of shear viscosity over 
entropy density for RHIC energies is fairly unconstrained. While the general trend of viscous 
corrections to ideal hydrodynamics has been studied by Teaney 



32] 



a dynamical implementation 
of viscous hydrodynamics and comparison to experimental data is still lacking. This is partly 
due to the fact that the "simplest" form of viscous hydrodynamics, the relativistic Navier-Stokes 
equations, are be-riddled by acausality problems and instabilities [33(1 . Therefore, there has been 



recent interest in so-called second-order (Israel-Stewart [34] ) theories 35, 0, S3, MM, H EE 5 
which are, however, of more complicated structure than the Navier-Stokes equations. 
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Specifically, it seems that adapting existing numerical hydro dynamic solvers to treat Israel- 
Stewart theory for all but the simplest geometries is a non-trivial task. It might therefore be worth- 
wile to devise completely new algorithms that are more suitable (or at least simpler) than present 
hydrodynamic solvers. Along these lines, in this work we present a straightforward algorithm 
for solving Israel-Stewart viscous hydrodynamic equations for geometries that are longitudinally 
expanding, are space-time rapidity independent and have radial symmetry and thus should be well 
suited to describe viscous hydrodynamics of central collisions at RHIC and in the future at the 
Large Hadron Collider (LHC). 

Our work is organized as follows: in section [XT] we review the equations of causal viscous hydro- 
dynamics and present evidence that our numerical algorithm reproduces the ideal hydrodynamic 
behavior in the limit of a small ratio of viscosity over entropy, as it should. 

In section IHIl we present a more involved test that is based on measuring correlation functions 
of small fluctuations and comparing to analytic results. 

In section HVl our results for the temperature evolution and particle spectra in relativistic causal 
viscous hydrodynamics are presented and we give our conclusions in section IVl 



II. SETUP AND COMPARISON WITH IDEAL HYDRODYNAMICS 

The basic equations of causal viscous hydrodynamics that we choose to study are given by (3^ 

(e + p)DvP = V»p - A^V CT IP CT + W v Bu v , (1) 
De = -(e + p)V^ + ^n^(V^ M ), (2) 
TnA£A£L>n Q/3 + IP" = r/( W) - 2r u U a ^u u l , (3) 

where e,p are the energy density and pressure, respectively, is the flow four-velocity that obeys 
u^u^ = 1 and IP" is the shear tensor that fulfills u^IP" = = 11^ and characterizes the viscous 
deviations in the energy momentum tensor, 

T*» = (e + pyV - pg^ + IP". (4) 

Furthermore, rj and tjj are the shear viscosity coefficient and relaxation time that are related by 
-2- = ^ in weakly-coupled QCD [39(] and the remaining definitions are 

d^u" = d^u" + T^u a , D = u^, V = A^d u , 

{A„B U ) = A^B U + A U B^ - l^ v A a B a , (A M , B v ) = \ (A^B U + A U B^) , (5) 

where are the Christoffel symbols. As outlined in the introduction, we will be interested in 
systems which are rapidity-invariant and have radial symmetry, therefore we choose to work in co- 
moving and radial coordinates r, r, <j>, 77 with the relations r = \/t 2 — z 2 , r 2 = x 2 + y 2 , tan <ft = y/x 
and r\ = atanh(z/t). The only non- vanishing fluid velocity components are then u T and u r with 
the relation u T = y/1 + (u r ) 2 and neglecting gradients in (ft and rj we find for the above equations 

(e + p)Du T = (l - (n T ) 2 ) {d T p - d v Il u T ) - u T u r (d r p - d v Ii u T ) 
{e+p)Du r = -u T u r {d T p-d u TV;)-(l + {u r ) 2 ){d r p-d u Tl») 

De = -( e + p)6 + ^ (-W r {\ - v 2 ) 2 {V r u r ) -r 2 nJ(V<V^) -r 2 nj(VV)) 
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-d v u v T = v 2 d T u r r + vd r u r r + n; [d T v 2 + d r v + ^ + ^ + In^ 
d„n u r = v d T n r r + d r u r r + n; [ d T v + - + + -m 

\ T T J T 

mDnv + iq = -r?r 2 (W) 

T U DU r r + U r r = -r]{V r u r ) + 2u r r u (Il r T Du T + U r r Du r ) 



(V r u r 



-2d r u r - 2u r Du r + § (1 + (u r ) 2 ) 9, 



-2^ + 1< 

r 3 



T 2 (V«'') 



e = d T u T + d r u r + + 



(6) 



where v = u r /u T , IT^ = — 1> II£, II? 



-IK — (1 — v 2 )W r and here D = u T d T + u r <9 r . This system of 



equations has to be closed by providing an equation of state, e.g. e = e(p). 

Clearly, it is possible to use the relation u T = yjl + (u r )' 2 to eliminate either u T or u r from the 
above equations. Defining 7 = u T = (1 — v 2 )~ l l 2 one obtains 



7 4 (6 + p)-(i- A 2 )n; 



- 7 2 (5 r + vd T ) P + (d r + w9 T )n; - 7 4 «(e + P ) + 7 2 ^n; 



(e+ P )7 2 -n; 



u9 T u — vd r e — (e + p) 



2 „ If 

7 2 d r u + - + - 
r r 



o r v 5- 



7} 



(7) 



However, in the code we prefer to keep both it 7 " and w r , solving equations for them independently 
so that a non-trivial consistency check on the numerics is provided by monitoring the deviation of 
(n r ) 2 — (u r ) 2 from unity. 



A. Discretization 

We use discretized space-time and compute differentials of a function /(r, r) as finite differences, 

»/(r,r) = /(T ' r + o) - /(T ' r - a) , ^/(r,r) = /(T + ^; ) - /(T ' r) , (8) 
la ot 

where a, 5r are the spatial and temporal lattice spacings, respectively. The boundaries are taken 
care of by using one-sided derivatives. 

Provided a starting condition at time r = tq for the variables u T , u r , e, II^, nj? one can then 
integrate the set of equations © forward in time. The virtue of this approach is that one imme- 
diately obtains the results for the fluid velocities etc. rather than having to perform the "usual" 
hydrodynamic algorithm (transforming to the calculational frame, integrating equations, trans- 
forming back). The drawback is that in Eqs. ©, time derivatives of the above variables are still 
coupled (e.g. the first equation of Eqs. © contains both d T u T and d T p). However, since all time 
derivatives enter only linearly this can be rectified by making use of a linear equation solver so 
that e.g. d T u T = f(r,r) which can be directly integrated using the above discretization. 

In practice, this works as follows: from the first equation in Eqs.([6j) we pick out the coefficients 
of the time derivatives d T u T , d T u r , d T p and label them as aoo, aoi> «02, respectively. The remaining 
part of the equation, which contains no time derivative, is called 6q- Thus, it becomes of the form 



ao d T u T + a m d T u r + a 02 d T p = b . 



0) 
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Note that in order to obtain this form we have expanded out the derivatives d^H", d u Ilf by using 
the relevant equations in ([6]). A similar procedure for the second and third equation of © leaves 
us with 

aiod T u T + a\id T u r + a\2d T p = b\ 

a2od T u T + a2id T u r + a 2 2d T p = b 2 . (10) 

With 5j = d T (u T ,u r ,p), these three equations may be written in matrix form as aijSj = bi, which 
has a solution unless detajj = 0. Numerically, this matrix equation is readily solved using a 
standard linear-equations solver, so 5j is known explicitly and may be used to finally compute 
d T in and d T W r . This completes the setup of our algorithm 1 . 



B. Testing the Code Ideal Hydrodynamics 



As a first test, we run our numerical code for a very small value of viscosity, r]/s = 10 and 
compare our results to ideal hydrodynamics. Our problem of choice is to start with a configuration 
for the energy density 



e(r, r ) 



e 



1 + exp [(r-Rj/o-y 



(11) 



where R = 6.4 fm can be thought of as the "radius" of a nucleus, and eo is such that the temperature 
(assuming an ideal gluon gas) is To = 0.2 GeV at r = 0. The parameter a is in principle arbitrary, 
but we choose it to be a = 0.02 fm in order to have a very steep fall of the energy density near 
r ~ R. Choosing the ideal equation of state e(p) = 3p for which the speed of sound squared (? s = |, 
we can then compare the time evolution of temperature and velocity to the analytic solution 

c?[l+(l-c s D JJ (r,r-ro))- 1 ]/2 



TRaym^, t) 

a R {r,t) 
VR(r,t) 



T e 



-c 8 a R (r,T-T ) ( T 



tanh 



a R (r,r - t ) + 



v R (r,T- r ) 



TO 



2 m V t+r 



+R 1- 



-R l+cs 





r~R+c s t 
t+c s (r-R) 

1 



OO 

r < R — c s t 
R-c s t<r<R + t 
r > R + t 



1 - c s v R (r, t 
r < R — c s t 
R-c s t<r<R + t 
r > R + t 



In 



TO 



TO 



(12) 



from Baym et al. [46J. Results are shown in FigQ3 where it can be seen that - within the errors of 
the approximate analytic solution - the numerical solution agrees with the results Eq.(|12p. 



III. FLUCTUATIONS AND LINEARIZED VISCOUS HYDRODYNAMICS 



Motivated by cosmolo gy where one can actually observe correlations of density fluctuations in 
the early universe 43|, 44, 45], we study radial fluctuations of the energy density e, the flow velocity 



1 A version of the code, written in reasonably well documented C, may be obtained upon request from 
paulrom@physik.uni-bielefeld.de 
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FIG. 1: Comparison between numerical results and the analytic approximation Eq. (fT2]) (full and dashed 
lines, respectively) for the evolution of the temperature (left figure) and the velocity v — ^ (right figure). 
The slight disagreement between numerical and analytical results is due to the approximations involved in 
the analytic solution and has the same sign and size as in the original work 46] . 



v and the shear tensor li^ around a background solution eo,u T ,H lu such that 

e(T,r) = e (T)+5e(T,r), v(r,r) = 6v(r,r), IP" = n^(r) + <5II^(t, r) 



where the background solution obeys the equations [3a, |47|] 
d T e Q -- 



T + T ii ^,o CtU^q- T ^r,,o+ g T 



(13) 



(14) 



and we remind that — = ^a. 

rn 3 

In what follows, we will assume that the fluctuations around the background are small so we 
keep only terms linear in the perturbations ("linearized hydrodynamics", c.f. [icj . |49|). Note, 
however, that we keep the full non-trivial time dependence of the background, which to the best 
of our knowledge has not been done before in the context of heavy-ion collisions even in the case 
of ideal hydrodynamics. 

To slightly simplify the discussion, we want to assume in this section that rn is constant with 
respect to time (which consequently requires a time-dependent ratio of r//s), whereas in other 
sections of this work tjj will be time-dependent. To linear order in the perturbations one is thus 
left with the set of coupled partial differential equations 



(e + Po - IE O )0 T + c 2 s d T e ~[d T + -) IE + -II? 



</.0 



Sv + c 2 s d r 5e 



8v + 



9 P0 



4 



d r + - 

r 

-2d r + - 
r 





d r + - 

r 


sn; 


1 + c 2 " 

T 


5e 


8 r 


2 r 






9 r 


d T + 


4c L 

— -5e + 
9 r 


d T + 



1 
1 



51% 



<5ie 







0. (15) 



These can be further simplified by noting that for the initial condition ITq = and no radial flow 
one has n^ = 11^ and as a consequence of 11^ = the relation IT^q = — ^11^ holds for all r. 
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Usually one would do a space-like Fourier transform to get rid of the space-like derivatives. Due 
to our choice of coordinates, however, this is obviously not possible. However, upon introducing 



SU = ( d r + - ) 8IZ + -<SIE» 



we can achieve the same goal by doing a so-called Bessel-Fourier transform, 

Sv(r,r) = Jq 00 dnJi(Kr)5v{T, k), Se(r,r) = J* °° dKJo(Kr)6e(r, k), 
5U(T,r) = / °°d«:Ji(/«r)<Jn(r,/s), SU^{T,r) = J °° dKj (Kr)Slty(T, tt) 

where the property of the Bessel functions J n 

drrJ n (Kr)J n (K'r) = — — 



l„n 2 „ 1 3 



(e + Po + ^lo)d T + cj d T e + - T + - 



Sv — Kc Q 8e — (511 



eo + Po + ^n^o 



k5v + 



1 + c 



5e - -5U2 








4 r~ 8c L~ 
- poK§v ---8e + 

8 2r~ 4 C L~ 

-p « dv - -K— de + 



d T + — 

m. 

d T + — 
m. 



5W n = 



<5n = o. 



(16) 



(17) 



(18) 



can be used to invert the above transformations. Using an ideal equation of state, po = c 2 eo and 
the Eqs. (|14p to remove explicit time derivatives on eo and 11^ we thus find 



(19) 



A. Sonic peaks in Ideal Hydrodynamics 

Upon first taking the limit m — > and then setting Big" as well as 5H^ 1 ' to zero we recover the 



equations for fluctuations in ideal hydrodynamics, which together with eo esc r 1 Cs require 



a 



2 _ 5 

r t 



d T + ^ + c 2 s k 2 



6v(t, k) = 



(20) 



(and a similar differential equation for 5e). As can be quickly verified, the solutions to the linearized 
ideal hydrodynamic equations then become 



5v (r, r) 
e>e(r, r) 



£teJi(Kr)r (1+c ' )/2 A(k) J ( _i +c 2 )/2 (kc s t) + B(k) Y(_ 1+c 2) /2 (kc s t) 



l + ci 



dKJo(Kr)T^ 1+c ^/ 2 \a{k) J {1+c 2 )/2 {kc s t) + B{k)Y {1+c 2 s)/2 {kc s t)] , (21) 



where J and Y are both Bessel functions of the first kind and A, B are constants of integration. 

As initial conditions at the starting time r = to we choose for simplicity 5v(jo,r) = and 
random noise for e)e(ro,r) with a correlation function 2 



e 2 < 5e(r ,r)5e(T ,r') >= A 2 — — — . 



(22) 



Let f^\r) be the i th configuration of an observable /. The correlation function < f(r)f(r') > is then defined as 
</(r)/(r')>=W^.»^E^i/ C0 W/ W (O- 
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which is the polar-coordinate equivalent of a Gaussian distribution with standard deviation A. 
This initial condition has the advantage that it implies 



< A(k)A(k') >= kS(k-k') 



such that 



A 2 c 



2^2 



^(l+ C 2)/2(«C s r ) - Y( 1+c 2y 2 (KC s T ) 



J(-l +c 2)/2{ K CsTo) 
Y (-l+cf)/2i Kc sTo) 



e 2 (r) < 5e(T,r)5e(T,r') >- 



dn k J (rer) J (Kr')f(K, r, r ) 



(23) 
(24) 



where 



f(K,T,T Q ) 



TO J 



A 2 



^(i+ c 2)/2( K c s r)y" ( _ 1+c 2) /2 (Kc s r ) - ■5 / "(i +c 2)/2(kc s t)J { _ 1+c 2) /2 (kc s to) 



2 • 



^(l+ C 2)/2( K C s To)Y"(_ 1+c 2)/ 2 (KC s r ) - y(l +c 2)/2(KC s To)J(_ 1+c 2)/ 2 (KC s To) 

Despite its ugly appearance, this is a nice result since for large k we find 

A 2 cos 2 (kc s (t — tq)), (26) 



(25) 



/(k,t,t ) -> — 

which are just the sonic peaks that one can also derive in cosmology. 

This result can serve as a stringent test on the numerical algorithm used to solve the hydro- 
dynamic equations as the position of the maxima and minima of / as a function of k are very 
sensitive to the speed of sound. In what follows, we thus initialize our numerical algorithm with 
precisely the same initial conditions as discussed above and then measure the correlation function 
e o~ ( r ) < Se(r,r) 5e(r,r') > to extract /(«,t,to) by using Eq. (fl~8j) to integrate out 3 both r,r', 



r dr r' dr' Jo(Kr)Jo(K'r')eQ (t) < 5e{r,r) 5e(r,r') > 



-f(K,T,T ). (27) 



To maximize the signal, we pick k = k' which is regular on the lattice we use to solve the hydro- 
dynamic equations. 

In FigJ2]we show the result for 4 /(ft, r, To) obtained on a lattice with lattice spacing a = 0.25 
GeV -1 and N = 8192 sites and rj/s = 10~ 4 , ensemble-averaged over 100 configurations and coarse 
grained in k, for three different times r. Up to three sonic peaks can be nicely distinguished and 
the comparison with the analytic result Eq. (|25p indicates that our code indeed accurately solves 
the ideal hydrodynamic equations with the "correct" speed of sound. There is, however, a slight 
discrepancy between the numerical measured correlation function and its analytic result at small 
k and later times: presumably, this is due to the fact that on the lattice, only a finite number 
of momenta can be simulated and thus the inversion formula Eq. (|27p holds only approximately. 
Indeed, the second part of Figl2]shows that this discrepancy can be systematically reduced by going 
to larger lattice volumes. Since the solution of the viscous hydrodynamic equations themselves do 
not depend on relations such as Eq. (|27l) . this discrepancy should not be mistaken as a failure of 
the algorithm to correctly treat low momentum modes. 



3 Since we solve the hydrodynamic equations on a lattice, in practice we do the integrations by summing over all 
lattice sites. Furthermore, momenta also discrete, where a is the spatial lattice spacing, N the number 
of lattice sites and k is a positive integer smaller than N/2. 

4 Note that the lattice dispersion relation k = a -1 sin (irk/N) has been used to convert to continuum values. 
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FIG. 2: Left figure: The correlation function /(rt,r, to) from solving hydrodynamic equations on a lattice 
(symbols, see text for details) for to = 1 fm/c and different t compared to the analytic result (full lines). 
As one can see, the agreement between the analytic result and the measured correlation function is very 
good in general (a slight difference e.g. in the speed of sound would be clearly visible by a shift in the 
minima/maxima of /). Note that at later times, a discrepancy at low momenta k develops. This is probably 
a lattice artifact since increasing the simulated volume reduces the discrepancy (right figure). 



B. Sonic Peaks in Viscous Hydrodynamics 

Treating the set of equations (|19p in their full generality we were unable to find analytic solutions 
like those obtained in the previous subsection. However, since together with Eq. (|14p these are just a 
set of six coupled ordinary differential equations they readily lend themselves to numerical solutions, 
which we nevertheless want to refer to as "analytic" since they are completely independent of our 
numerical algorithm to solve the hydrodynamic equations. 

With the same initial conditions as in the previous subsection one therefore obtains an "analytic" 
solution of the correlation function /(k, t,to), with both expansion and viscosity included. In 
FigGJ this "analytic" solution is again compared to the correlation function obtained by solving 
the hydrodynamic equations on a lattice (for r//s = 0.1 and rj/s = 0.3, respectively). Similarly 
to the case of vanishing viscosity, we find that there is very good general agreement between the 
measured (ensemble averaged and coarse-grained) correlation function and the "analytic" result, 
except for later times and small momenta k, where lattice artifacts seem to be accumulating. 

Since also in this case the position and width of the sonic peaks are very sensitive to the value of 
rj/s and the speed of sound, we argue that the good agreement between measured and "analytic" 
correlation functions is a strong indication that our numerical code is indeed correctly solving the 
second order viscous hydrodynamic equations. 

Finally, we want to point out that fluctuation measurements may also help to constrain the 
value of r//s from RHIC data, as has been recently suggested [50| . 



IV. CAUSAL VISCOUS HYDRODYNAMICS WITH TRANSVERSE FLOW 

Let us now study the effects of viscosity on quantities of interest for heavy-ion collisions. For 
simplicity, we will assume that the radial energy density profile is given by Eq. (|Llj) . where we take 



Rq = 6.4 fm and a = 0.54 fm which has been used before for ideal hydrodynamic calculations [461 ] . 
The constant eo is chosen such that we have an initial temperature To at r = 0. In accordance 
with ideal hydrodynamic studies in their simplest form we assume that at the time when we start 



9 




FIG. 3: Comparison for /(/t, t,tq) from solving hydrodynamic equations on a lattice (symbols, rj/s = 0.1 
(left) and rj/s — 0.3 (right), respectively) with 4096 sites and lattice spacing a — 0.25 GeV -1 and "analytic" 
solution (full lines) of Eqs.JTH). 



applying our hydrodynamic description (r = tq), the system does not have any transverse flow 
already, so v = 0. Furthermore, since this is still an exploratory study, we pick a simplistic equation 
of state, p = c 2 s e with constant speed of sound c 2 s = 1/3. 

This set of initial conditions would be sufficient to determine the subsequent dynamics fully in 
the case of ideal hydrodynamics. Including the effects of viscosity requires that we pick a specific 
value of the ratio rj/s and also provide initial conditions for the two independent components of 
IF" 1 ' at r = to- Maybe the simplest choice would be to assume - like in ideal hydrodynamics - 
that the system for some reason happens to be in equilibrium at r = to, such that "accidentally" 
IP 1 ' = 0. This choice probably highlights best the difference of viscous hydrodynamics to ideal 
hydrodynamics, since one starts from the same initial condition, so we will use it as our initial 
condition in the following. 

However, there are other "sensible" choices of initial conditions that might be more relevant 
for real heavy-ion collisions in the future. E.g. within the Color-Glass-Condensate model in its 



simplest form (the McLerran-Venugopalan model 51f] ) , the system does not have any longitudinal 



dynamics, so after times r > Q s 1 where Q s is the saturation scale, the system essentially has zero 



longitudinal pressure 52|, [53j. In the local rest frame u r = = vP = 0, so with c 2 s = 1/3, Eq.(jH) 
would imply 

U^=p, K = ~~ (28) 

Finally, going beyond the McLerran-Venugopalan model to include so called next-to-leading 



order corrections of gluon production [54|, [55|] one has to take into account the effect of rapidity 
fluctuations and full three-dimensional gauge field dynamics. This has recently been shown to 
trigger instabilities [3(J, leading to the generation of non-zero longitudinal pressure [56| at r > Qj 1 . 
Pending the result using the correct rapidity fluctuation spectrum [HtJ , the initial condition for IT"' 
is expected to lie somewhere in-between the two cases discussed above. 



A. Temperature Profile in Viscous Hydrodynamics 

Choosing the initial condition = at r = To, we can investigate the changes of the 
temperature profile from the ideal hydrodynamic behavior due to dynamical viscous effects. In 
FigHJwe show the temperature as a function of the radius for Tq = 0.36 GeV at to = 1 fm/c, but 
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FIG. 4: Temperature profile for calculations with different rj/s (dashed, dotted and solid lines, respectively) 
for three different times (see text for details). As expected, for larger values of rj/s, differences to ideal hydro- 
dynamics are biggest and viscous hydrodynamics initially cools slower than ideal hydrodynamics. However, 
note that in certain regions and at later times, viscous hydrodynamics turns out to give temperatures slightly 
smaller than the corresponding ideal hydrodynamic calculation. 

different values of rj/s, calculated on a lattice with 512 sites and a = 0.25 GeV -1 lattice spacing. 
Choosing the temporal time step as 5t = 0.005a we find that the violation \J\{u T ) 2 — (u r ) 2 — 1|, 
summed over all lattice sites and divided by the number of sites always stays smaller than one 
percent, providing yet another check on the numerics. Finally, we have checked that choosing a = 1 
or 0.5 GeV -1 does not result in any noticeable deviations of our calculated temperature profile, 
thus we have some confidence that our results are not strongly affected by numerical artifacts. 

For early times, the behavior shown in FigJH shows that - as in the case of neglecting transverse 
flow [i^, [33] - in viscous hydrodynamics the temperature decreases slower than in ideal hydrody- 
namics. However, at late times this behavior is seemingly counteracted by viscous radial dynamics: 
at very small values of the radius, the viscous hydrodynamic equations result in slightly lower tem- 
peratures than in the ideal hydrodynamic case. Note that such behavior has been also found in 



The success of the hydrodynamic picture in the context of heavy-ion collision builds upon 
the ability to fit the particle spectra observed in these collisions. While now methods of how to 
convert hydrodynamic quantities such as energy density and fluid velocity into particle spectra 
have reached some sophistication, the main building block still seems to involve the Cooper-Frye 
freeze-out prescription [60(] in some form or the other, which states that the spectrum of particles 
with energy E and momentum p is given by 



Sinn. 



B. Particle Spectra in Viscous Hydrodynamics 




(29) 
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FIG. 5: Inverse slope parameter T s i ope for Tf = 0.135 GeV and tq = 1 fm/c. Two initial conditions for IP", 
corresponding to pressure isotropy (full line) and vanishing longitudinal pressure (dashed line) at r = To are 
shown. Choosing To = 0.36 GeV (left), the spectra become increasingly flatter when raising r]/s, while this 
effect can be compensated by lowering Tq (right, shown for rj/s = 0.16). 



where d is the degeneracy of the particles and is the velocity which comes out of the solution 
of the hydrodynamic equations. Here / is the distribution function which - including viscous 
corrections - can be written as 



32, 3 



' = A ( 1 + #tTri)' (30) 

where for simplicity we take the equilibrium distribution /o to be given by the Boltzmann- 
distribution 5 fo(x) = exp(— x). 

Furthermore, cffi' 1 is the normal vector on the freeze-out surface, parametrized as dS^ = 
^cosh 7], cos <p —^^ , sin <ft ~^~ ; sinh r^J rdrTj(r)d(j)dri in our choice of coordinates |37.]. Here T/(r) 
is the freeze-out time parametrized as a function of r or - put differently - the time at which the 
slab of matter at radius r has reached the freeze-out condition. 

For this exploratory study, we will apply the Cooper-Frye freeze-out prescription to convert 
the hydrodynamic variables at a single specific temperature (the freeze-out temperature Tf) into 
transverse momentum spectra for particles. This is what also has been used in early ideal hydrody- 



namic calculations [611 . |62| . Since we use a gluonic equation of state and do not include a realistic 
matching to hadronic degrees of freedom, we contend ourselves to study the effects of viscosity on 
the spectrum of gluons mostly. 

For the spectrum, in terms of particle transverse momentum p±, angle 4> p and rapidity y one 
thus finds 

/ drf(r)\ 
Pfj.dT,' 1 = f m_|_ cosh(?7 — y) —p_\_cos{4> — 4>p) — j j Tf(r)rdrd(j)dr], (31) 



with m± = \Jp\_ + wiq and mo the rest mass of the particle. 

Since in our calculations we are only including transverse flow we have 

Pfj,u^ = (m± cosh(r/ — y)u T — p± cos(^ — <p p )u r ) , (32) 



5 It is easy to change this to Bose-Einstein or Fermi-Dirac distributions, but for massive particles we have found 
the differences in resulting observables to be minimal. 
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FIG. 6: Mass dependence of viscous effects: shown are spectra for pions, kaons and protons for a freeze-out 
temperature of Tt = 0.135 GeV and for To = 0.36 GeV and IP" = at To = 1 fm/c. As can be seen, the 
higher the mass of the particle, the stronger do viscous effects change its low p± behavior. 



which allows us to integrate out both angles <f> and rapidities rj in Eq.([29|). Using 

d 

dx , 

r* dtj> 



drj cosh n 7yexp(— x cosh??) = 
/ — — cos n<f) exp (x cos (f>) = I n (x), 

JO TT 



2K (x 



where K n (x) and I n {x) are modified Bessel functions one finds for the particle spectrum 

d 3 N d 3 N d 3 8N 
E — — — = E^A + E 



d 3 p d 3 p 
with the equilibrium part taking the form 
d 3 N 2d 



E 



d 3 p (2tt) 2 



rdrTf(r) 



m±I (u r p±/T)K 1 (u T m ± /T) 



d 3 p ' 
dT f {r) 



(33) 



(34) 



dr 



P±h(u r p±./T)K (u T m ± /T) 



(35) 

where the integral over r runs from to the maximum freeze-out radius if Tf(r) is a single- valued 
function (else one has to introduce a different parametrization of the freeze-out surface). Noting 
that the Bessel K functions always have the argument u T m±/T (and similarly for the J's) we 
refrain from writing the argument in the following. Noting that 



p^PvU.^ = H T T 2vm±p± cosh(y — rj) cos(^ p — (f)- 
—v 2 m 2 L cosh 2 (y — rj) — v 2 p\ s\r?(<p p - 
+11^ -m\ sinh 2 (y - rj) +p± sm 2 (4> p 
we then find the for the dissipative corrections to the spectrum 



p\ cos(2(0p - <j))) 



(36) 



E- 



d 3 5N 



d 



d 3 p (2tt) 2 



rdr 



Tf(r) 



2T 2 (e+p) 



2p±m±vIL r r 



drt 



m ± (K + K 2 )h ~ p±-r L K 1 (I 2 + I ) 
dr 
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-pIK 



2m ± K 1 I 2 - P± -lKo(h + h: 



-v 2 m 2 ± U r r 



-vV±K 



dr 

^m ± (3K! + K a )I - p±^ l (K + K 2 )h 
2 dr 

- I 2 ) - pJ^lKvUh ~ h) 
dr 2 



+p\^ 
-m 2 , m 



m ± K 1 (I - I 2 ) - p J^lK Q \(h - h) 
dr 2 

m ± ±(K 3 - K X )I - P± ^(K 2 - K )h 
2 dr 



(37) 



where we remind that v = u r /u T . 

In FigEJ we show the inverse slope parameter T s i ope of gluons which we define by calculating 
the gluonic spectrum at Tj = 0.135 GeV and fitting it by 



E 



,d 3 N 
d?p 



1 



T 2 

slope 



exp [-p±/T s i ope ] 



for 0.2 < p ± < 1 GeV (c.f.[63|]). In FigG^a), the slope has been calculated for T = 0.36 GeV, 
r = 1 fm/c and for the two extreme cases H^ v {tq) = (full line) and zero initial longitudinal 
pressure Eq. (|28p (dashed line). As can be seen from this figure, increasing r]/s and leaving all other 
parameters unchanged leads to an increasing T s i ope ( "flatter spectra" ) for gluons, with no dramatic 
difference between the two different initial choices for IP* 7 . However, as has been anticipated from 
our earlier studies neglecting the effect of transverse flow [3^], one can compensate this effect by 
changing the effective initial conditions. This can be seen in Figj5jb), where we show the spectral 
slope for the same freeze-out temperature, but different initial temperatures To. 

It is also interesting to study how the presence of viscosity affects massive particles. To this 
end, hypothetical spectra of pions, kaons and protons for Tf = 0.135 GeV, To = 0.36 GeV and 
n^ ll/ ( T o) = at ro = 1 fm/c are shown in Figj6j These spectra cannot be directly interpreted as 
real particle spectra because a realistic matching to a hadronic equation of state and the effects 
from higher mass resonance decays 64j] are missing in this study 6 . Nevertheless, from Fig E] one 
can glean that the more massive a particle is, the more viscosity affects its spectrum, in particular 
at low p±. Indeed, this can be traced back to Eqs. (135|37p which in the limit of vanishing p± and 
neglecting radial dynamics (v = 0) predict negative d?N/d?p for large tuq/T, more specifically for 



mo 2(e + p) - 
T > UV 



—IT 1 

8 1L V 



(38) 



'/ 



Thus it seems that - whenever 11^/ (e + p) becomes non-negligible - viscous corrections 5N to the 
spectrum of high-mass particles become very large, e.g. more than 100 percent at low p±. While 
it is unclear at which value of rj/s this starts to be a problem in practice, it nevertheless serves 
as an indication that the assumption of small deviations from equilibrium [3^] is breaking down. 
Consequently, the reliability of the tool we have used to probe the system dynamics, namely viscous 
hydrodynamics, becomes questionable. Thus, for r]/s larger than a critical value, one probably has 
to use a different model than hydrodynamics to correctly calculate observables that are to be 
compared to experiments. 



See however [65[ for a comparison to experimental data. 
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V. CONCLUSIONS 

We have studied the effect of shear viscosity in a hydrodynamic description of central heavy-ion 
collisions. We presented a simple algorithm to solve the relevant equations numerically and have 
successfully carried out several tests on this algorithm. These tests are not specific to our algorithm, 
but can in general be used to test any algorithm for solving relativistic viscous hydrodynamics. 
Assuming an ideal equation of state e = 3p for simplicity we calculated the time evolution of the 
temperature profile of a central heavy-ion collision, finding that while viscous hydrodynamics in 
general cools slower, certain regions at later times may cool faster than in a corresponding ideal 
hydrodynamic calculation. 

We also calculated the effect of viscosity on the slope of gluon spectra, finding that for small 
values of rj/s, changes can largely be compensated by lowering the temperature at which the hy- 
drodynamic evolution is started. For massive particles we find that viscosity changes the spectrum 
the more the higher the mass of the particle under consideration. We give arguments that for a 
sufficiently large value of 77/s, corrections which in the derivation of the viscous hydrodynamic equa- 
tions had been assumed to be small actually become large, thus signaling the possible breakdown 
of any hydrodynamic description of the system. 

Even though our simplifying assumptions (ideal equation of state, no feed-down correction, 
only radial flow) leave ample room for improvement, we hope that our study provides the basis for 
coming viscous hydrodynamic algorithms as well as fits to experimental data. 
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